/**
 * 
 */
package org.estela.stats;

import org.estela.util.Math2;


/**
 * Simple kernel density estimator. 
 * Uses one gaussian kernel per observed data value.
 * 
 * @author Antonio Kantek
 */
public class KernelEstimator implements ProbabilityEstimator {

	/** Maximum percentage error permitted in probability calculations */
	private static double MAX_ERROR = 0.01;

	/** Values seen */
	private double[] _values;

	/** Associated weights */
	private double[] _weights;

	/** Number of values stored in _weights and _values so far */
	private int _numValues;

	/** Sum of the weights */
	private double _sumOfWeights;

	/** Values standard deviation */
	private double _stddev;

	/** Precision of data values */
	private double _precision;

	/** Whether we can optimise the kernel summation */
	private boolean _allWeightsOne;

	/**
	 * 
	 */
	public KernelEstimator(double precision) {
		assert precision > 0 : "precision cannot be zero";
		_values = new double[50];
		_weights = new double[50];
		_numValues = 0;

		_sumOfWeights = 0;
		_allWeightsOne = true;
		_precision = precision;

		if (precision < Math2.SMALL_CONSTANT) {
			_precision = Math2.SMALL_CONSTANT;
		}
		_stddev = _precision / (2 * 3);
	}

	private int findNearestValue(double key) {

		int low = 0;
		int high = _numValues;
		int middle = 0;
		while (low < high) {
			middle = (low + high) / 2;
			double current = _values[middle];
			if (current == key) {
				return middle;
			}
			if (current > key) {
				high = middle;
			} else if (current < key) {
				low = middle + 1;
			}
		}
		return low;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see com.dme.stats.Estimator#addValue(double, double)
	 */
	public void addValue(double data, double weight) {
		if (weight == 0) {
			return;
		}
		data = Math.rint(data / _precision) * _precision;
		;
		int insertIndex = findNearestValue(data);
		if ((_numValues <= insertIndex) || (_values[insertIndex] != data)) {
			if (_numValues < _values.length) {
				int left = _numValues - insertIndex;
				System.arraycopy(_values, insertIndex, _values,
						insertIndex + 1, left);
				System.arraycopy(_weights, insertIndex, _weights,
						insertIndex + 1, left);

				_values[insertIndex] = data;
				_weights[insertIndex] = weight;
				_numValues++;
			} else {
				double[] newValues = new double[_values.length * 2];
				double[] newWeights = new double[_values.length * 2];
				int left = _numValues - insertIndex;
				System.arraycopy(_values, 0, newValues, 0, insertIndex);
				System.arraycopy(_weights, 0, newWeights, 0, insertIndex);
				newValues[insertIndex] = data;
				newWeights[insertIndex] = weight;
				System.arraycopy(_values, insertIndex, newValues,
						insertIndex + 1, left);
				System.arraycopy(_weights, insertIndex, newWeights,
						insertIndex + 1, left);
				_numValues++;
				_values = newValues;
				_weights = newWeights;
			}
			if (weight != 1) {
				_allWeightsOne = false;
			}
		} else {
			_weights[insertIndex] += weight;
			_allWeightsOne = false;
		}
		_sumOfWeights += weight;
		double range = _values[_numValues - 1] - _values[0];
		if (range > 0) {
			_stddev = Math.max(range / Math.sqrt(_sumOfWeights), _precision
					/ (2 * 3)); // allow at most 3 sds within one interval
		}
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see com.dme.stats.Estimator#getProbability(double)
	 */
	public double getProbability(double data) {
		double delta = 0, sum = 0, currentProb = 0;
		double zLower = 0, zUpper = 0;
		if (_numValues == 0) {
			zLower = (data - (_precision / 2)) / _stddev;
			zUpper = (data + (_precision / 2)) / _stddev;
			return (Probability.normalProbability(zUpper) - Probability
					.normalProbability(zLower));
		}
		double weightSum = 0;
		int start = findNearestValue(data);
		for (int i = start; i < _numValues; i++) {
			delta = _values[i] - data;
			zLower = (delta - (_precision / 2)) / _stddev;
			zUpper = (delta + (_precision / 2)) / _stddev;
			currentProb = (Probability.normalProbability(zUpper) - Probability
					.normalProbability(zLower));
			sum += currentProb * _weights[i];
			/*
			 * System.out.print("zL" + (i + 1) + ": " + zLower + " ");
			 * System.out.print("zU" + (i + 1) + ": " + zUpper + " ");
			 * System.out.print("P" + (i + 1) + ": " + currentProb + " ");
			 * System.out.println("total: " + (currentProb * m_Weights[i]) + "
			 * ");
			 */
			weightSum += _weights[i];
			if (currentProb * (_sumOfWeights - weightSum) < sum * MAX_ERROR) {
				break;
			}
		}
		for (int i = start - 1; i >= 0; i--) {
			delta = _values[i] - data;
			zLower = (delta - (_precision / 2)) / _stddev;
			zUpper = (delta + (_precision / 2)) / _stddev;
			currentProb = (Probability.normalProbability(zUpper) - Probability
					.normalProbability(zLower));
			sum += currentProb * _weights[i];
			weightSum += _weights[i];
			if (currentProb * (_sumOfWeights - weightSum) < sum * MAX_ERROR) {
				break;
			}
		}
		return sum / _sumOfWeights;
	}

	public String toString() {

		String result = _numValues + " Normal Kernels. \nStandardDev = "
				+ _stddev + " Precision = " + _precision;
		if (_numValues == 0) {
			result += "  \nMean = 0";
		} else {
			result += "  \nMeans =";
			for (int i = 0; i < _numValues; i++) {
				result += " " + _values[i];
			}
			if (!_allWeightsOne) {
				result += "\nWeights = ";
				for (int i = 0; i < _numValues; i++) {
					result += " " + _weights[i];
				}
			}
		}
		return result + "\n";
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		try {
			KernelEstimator newEst = new KernelEstimator(0.01);
			newEst.addValue(1, 1);
			newEst.addValue(2, 1);
			newEst.addValue(3, 1);
			newEst.addValue(1, 5);
			System.out.println(newEst);

			double start = 1;
			double finish = 2;
			for (double current = start; current < finish; current += (finish - start) / 50) {
				System.out.println("Data: " + current + " "
						+ newEst.getProbability(current));
			}
		} catch (Exception e) {
			System.out.println(e.getMessage());
		}

	}

}
